rm(list=ls())
set.seed(1987)
install.packages('foreign')
library(foreign) 

# Insert path to main folder
#setwd("C:/replication")
setwd("//osiride-fs/m030221/private/ricerca/maj status/replication")


### In this R file are collected the codes for generating samples used in the simulation
### section of the paper (Section 3). 
### In particular, we have two possible simulations: the main one presented in the paper 
### (and detailed in supplement D),
### and the alternative simulation (of supplement E) 


# MAIN SIMULATION (SECTION 3) ---------------------------------------------
###  Parameter settings
n <- 601 # Number of districts
anni <- 100 # Number of years
LEAD <- runif(anni) # Draw randomly the identity of the winning party 

### Generation of the samples
dat <- NULL
yr <- rnorm(1,0,0.05)
for(t in 1:anni){
  Dt <- matrix(ncol=16,nrow=n)
  
  # Draw the identity of the winning party 
  lead <- LEAD[t]
  
  # Generate the hyperparameter of the Beta distribution needed later 
  if(lead>0.5){# if democrats hold control
    theta<-runif(1,5.1,5.5) 
  }else{# if republicans hold control
    theta<-runif(1,4.5,4.9) 
  }
  
  # Generate the vote shares in each district
  X <- rbeta(n,theta,10-theta)
  
  # Create variables 
  Dt[,1] <- rep(t,n) # year
  Dt[,2] <- 1:n # district id
  Dt[,3] <- X # vote share 
  Dt[,4] <- X-0.5 # winning gap (running variable)
  Dt[,5] <- ifelse(X>0.5,1,0) # D: dummy variable 1 if democrats win in that district-year, 0 otherwise
  Dt[,6] <- ifelse(mean(ifelse(X>0.5,1,0))>0.5,1,0) # dummy variable 1 if democrats have control of the assembly in that year
  Dt[,7] <- Dt[,4]*Dt[,4] # winning gap^2
  Dt[,8] <- Dt[,4]*Dt[,5] # interaction term between running variable and D
  Dt[,9] <- Dt[,8]*Dt[,8] # interaction term between running variable and D squared
  Dt[,10] <- ifelse(Dt[,6]==Dt[,5],1,0) # M: dummy 1 if distict belongs to the majority in that year
  Dt[,11] <- Dt[,5]*Dt[,10] # interaction between D and M
  # different DGPs of Y below
  Dt[,12] <- log(X+1) 
  Dt[,13] <- 20*X*(X*X-X+0.1)+0.5
  Dt[,14] <- 0.5+0.3*Dt[,5]+0.3*Dt[,10]+0.5*Dt[,6]+yr+rnorm(n,0,0.03)+Dt[,13]
  Dt[,15] <- 0.5+0.3*Dt[,5]+rnorm(n,0,0.03)+Dt[,13]
  Dt[,16] <- 0.5+0.3*Dt[,5]+0.3*Dt[,10]+rnorm(n,0,0.03)+Dt[,13]
  
  dat<-rbind(dat,Dt)
  
}

colnames(dat)=c("t","id","ds","dgap","d",'dm','dgap2','rdgap','rdgap2','m',
                'inter_d_m','log_ds','poly','y_full','y_d','y_d_m')

### Save data in .dta
write.dta(as.data.frame(dat), paste0("data/figure1/simulation_main.dta"))




# ALTERNATIVE SIMULATION (APPENDIX) ---------------------------------------
rm(list=ls())
set.seed(2022)

###  Parameter settings
n <- 601 # Number of districts
anni <- 100 # Number of years
LEAD <- runif(anni) # Draw randomly the identity of the winning party 
z <- 30 # Average number of competitive districts (choose among 30, 60, 90 which correspond to approx 5%, 10% and 15%)

### Generation of the samples
dat <- NULL
yr <- rnorm(1,0,0.05)
  
for(t in 1:anni){
    Dt <- matrix(ncol=16,nrow=n)
    lead <- LEAD[t]
    
    # Generate the hyperparameter of the Beta distribution needed later 
    if(lead>0.5){# if democrats hold control
      theta<-runif(1,0.51,0.55) 
    }else{# if republicans hold control
      theta<-runif(1,0.45,0.49) 
    }
    
    # Generate competitive districts and their vote share
    if(z==30){n1 <- 51}else if(z==60){n1<-101}else if(z==90){n1=151}
    theta <- 100*theta
    fat <- 100
    X1 <- rbeta(n1,theta,fat-theta)
    
    # Generate Republican districts and their vote share
    nrep <- (n-n1)/2
    theta <- 150
    fat <- 400
    X2 <- rbeta(nrep,theta,fat-theta)
    
    # Generate Democratic districts and their vote share
    ndem <- nrep 
    theta <- 250
    fat <- 400
    X3 <- rbeta(ndem,theta,fat-theta)
    
    X<-c(X1,X2,X3)
    
    # Create variables 
    Dt[,1] <- rep(t,n) # year
    Dt[,2] <- 1:n # district id
    Dt[,3] <- X # vote share 
    Dt[,4] <- X-0.5 # winning gap (running variable)
    Dt[,5] <- ifelse(X>0.5,1,0) # D: dummy variable 1 if democrats win in that district-year, 0 otherwise
    Dt[,6] <- ifelse(mean(ifelse(X>0.5,1,0))>0.5,1,0) # dummy variable 1 if democrats have control of the assembly in that year
    Dt[,7] <- Dt[,4]*Dt[,4] # winning gap^2
    Dt[,8] <- Dt[,4]*Dt[,5] # interaction term between running variable and D
    Dt[,9] <- Dt[,8]*Dt[,8] # interaction term between running variable and D squared
    Dt[,10] <- ifelse(Dt[,6]==Dt[,5],1,0) # M: dummy 1 if distict belongs to the majority in that year
    Dt[,11] <- Dt[,5]*Dt[,10] # interaction between D and M
    # different DGPs of Y below
    Dt[,12] <- log(X+1) 
    Dt[,13] <- 20*X*(X*X-X+0.1)+0.5
    Dt[,14] <- 0.5+0.3*Dt[,5]+0.3*Dt[,10]+0.5*Dt[,6]+yr+rnorm(n,0,0.03)+Dt[,13]
    Dt[,15] <- 0.5+0.3*Dt[,5]+rnorm(n,0,0.03)+Dt[,13]
    Dt[,16] <- 0.5+0.3*Dt[,5]+0.3*Dt[,10]+rnorm(n,0,0.03)+Dt[,13]
    
    dat<-rbind(dat,Dt)
    plot(density(X))
    
  }
  
colnames(dat)=c("t","id","ds","dgap","d",'dm','dgap2','rdgap','rdgap2','m',
                  'inter_d_m','log_ds','poly','y_full','y_d','y_d_m')

### Save data in .dta
write.dta(as.data.frame(dat), paste0("data/figure1/simulation_alternative.dta"))
  




